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I. INTRODUCTION 

State of the art lattice calculations of hadronic ma- 
trix elements have produced very accurate spectroscopic 
information. Examples of the accuracy reached in 
quenched lattice QCD are the calculation of the masses of 
low lying hadrons [1] and glueballs [2] . However progress 
in determining hadron wave functions and quark distri- 
butions has not been so rapid. The initial calculations of 
Bethe-Salpeter amplitudes were carried out on a rather 
small lattice in the mid 80's [3] for the pion and the rho. 
Further progress came in the early 90's when the Bcthc- 
Salpetcr amplitudes were calculated on a larger lattice 
for the pion and the rho [4,5] as well as for the nucleon 
and the A [4]. However the results were of limited in- 
terest because of their manifest dependence on the gauge 
chosen. A different approach to explore hadronic struc- 
ture was pursued by the authors of refs. [6]; instead of 
fixing a gauge or a path for the gluons, they considered 
correlation functions of quark densities which, being ex- 
pectation values of local operators, are gauge invariant. 
This is the approach we have adopted in this work. 

Our main motivation for studying density correlators is 
that they reduce, in the non-relativistic limit, to the wave 
function squared and thus they provide detailed, gauge- 
invariant information on hadron structure. The shape of 
hadrons is one such important quantity that can be di- 
rectly studied. The issue whether the nucleon is deformed 
from a spherical shape was raised twenty years ago [7] and 
is still unsettled. Because the spectroscopic quadrupole 
moment of a spin one-half particle vanishes, in experi- 
mental studies one searches for quadrupole strength in 
the 7*iV — -> A transition. Spin-parity selection rules 
allow a magnetic dipole, Ml, an electric quadrupole, E2, 
or a Coulomb quadrupole, C2, amplitude. If both the 
nucleon and the A are spherical then the electric and 



Coulomb quadrupole amplitudes are expected to be zero. 
Although Ml is indeed the dominant amplitude there is 
mounting experimental evidence that E2 and C2 are non- 
zero. The physical origin of a non-zero E2 and C2 am- 
plitude is attributed to different mechanisms in the var- 
ious models. In quark models the deformation is due to 
the colour-magnetic tensor force [7]. In "cloudy" baryon 
models it is due to meson exchange currents [8]. 

A recent experimental search at q 2 = 0.126 GeV 2 has 
yielded an electric quadrupole to magnetic dipole ampli- 
tude ratio [9]: 



Rem = WSmi = (-2.1 ± 0.2 ± 2.0)% 



(1) 



The larger error, an order of magnitude larger than the 
statistical error, is due to the model dependence in the 
extraction of this ratio from the experimental data. The 
experimental determination of Rem is complicated by 
the presence of nonresonant processes coherent with the 
resonant excitation of A(1232). 

The ratio Rem can be evaluated within lattice QCD 
without any model assumptions by computing the transi- 
tion matrix element N to A. An early lattice calculation 
of this transition matrix element provided an estimate of 
#bm = (-3±8)% [10]. 

Here we consider an alternative route to understanding 
the issue of deformation, via the direct study of hadron 
wave functions. We compute density-density correlators 
for mesons and baryons and three-density correlators for 
baryons and look for asymmetries when these are pro- 
jected along the spin axis or perpendicular to it. The 
observables that we use are described in section II. In sec- 
tion III we give the relations for the deformation among 
the states of different spin projections. Our quenched 
lattice results are presented in section IV. They clearly 
give support to a deformed rho. The nucleon deforma- 
tion averages to zero in agreement with the fact that its 



1 



spectroscopic quadrupole moment is zero. An analysis of 
the intrinsic nucleon deformation requires determining 
the body-fixed coordinates by diagonalization of the mo- 
ment of inertia tensor, which must be done configuration 
by configuration. This is too noisy to yield a statistically 
significant result. On the other hand the A has a non 
zero spectroscopic moment and any deformation should 
be detected via projection with respect to its spin axis. 
In quenched QCD we detect no significant deformation 
for the A within our present statistics. 

In addition to providing information about deforma- 
tion, baryon wave functions, calculated for the first time 
in a gauge invariant way for all values of the relative co- 
ordinates, can be used to study indirectly the potential 
among the three quarks. By performing fits to the three- 
density correlators, we find that a reasonable ansatz for 
the baryonic potential is provided by the sum of two-body 
potentials, called the A- Ansatz [11]. 

Comparison between quenched and full QCD results 
sheds light on the role of the pion cloud. It is known 
that quenching eliminates all or part of the pion cloud 
depending on the hadronic state. For the rho channel 
no intermediate backgoing quarks with the pion quan- 
tum numbers are present [12], so the deformation that 
we observe is not due to the pion cloud. We use the 
SESAM configurations [13] to investigate unquenching ef- 
fects. We find that the deformation in the rho increases. 
We detect a deformation in the A which remains small 
for the values of the dynamical quark mass considered. 
The unquenched lattice results are presented in section 
V. Section VI contains our conclusions. 



II. CORRELATION FUNCTIONS 

The wave function of a meson is usually defined in a 
given gauge g as the equal time Bethe-Salpeter amplitude 



$f(r) = J d 3 r' < 0\q fl (r')Tq f2 (r' +r)\M 



> 



(2) 



Summation over r' projects onto zero momentum. For 
the A + , the interpolating field that we take is 



where T is a Dirac matrix with the quantum numbers of 
the meson M of flavor /i, fa. <I>^ s (r) is the minimal Fock 
space state wave function, which is an approximation to 
the full wave function since other multiquark components 
are excluded. For the Bethe-Salpeter amplitude one must 
either fix the gauge or connect the quarks with gluons 
to form gauge-invariant (but path-dependent) quantities. 
Wave functions for baryons are defined in an analogous 
way to the meson wave functions but they involve two 
relative distances. E.g. the Bethe-Salpeter amplitude, 
<f>^ s (ri, r 2 ), for the proton is given by 

< 0\e abc u a s (r', t) (u b T (r' + r u t) C l5 d c (r' + r 2 , t)^j \B > 

(3) 



^abc 



« a (x) 2u bT (x) C 7 ^ c (x) 



+ d a (x)(V T (x) C 7 ,« c (x)) 



(4) 



where C = 7072 is the charge conjugation operator. Wc 
use r± = (71 =F ^72) /2 to create a hadron of a definite 
spin component J z . Explicitly for the first term in Eq. 4 
we take 



h/2 



^abc 

7f 

~abc 



Jl/2 : 7f 
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<(x)(2 U fcT (x) C 73 d c (x) 
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J - 1/2 : 7f 



°(x)(2u feT (x) C 73 rf c (x) 



+<(x)(2w bT (x) CT_d c (x) 



abc 



J-3/2 — — 



V3 



■5(x)( 2u bT {x) CF_d c (x 



(5) 



where Qi is the zth component of the spinor. The same 
construction is done for the second term of Eq. 4. With 
these combinations we recover, in the non-relativistic 
limit, the quark model wave functions of these states. 

The Bethe-Salpeter amplitudes in the Coulomb and in 
the Landau gauge were investigated in ref. [4]. Instead 
of fixing a gauge, a gauge-invariant wave function, which 
corresponds to the Bethe-Salpeter amplitude in the axial 
gauge, can be constructed by joining the quark and the 
antiquark with a thin string of glue [14]: 



My) = <o|g(o)r e *r ' dzA ^ q { y )\ M ) 



(6) 



Other variations to the thin string are to smear the glu- 
ons [5] or to evolve them so that they reach their ground 
state distribution. It was shown in ref. [14] that there 
is a much larger probability to find a quark-antiquark 
pair separated by, say, 1 fm when connected by a phys- 
ical adiabatic flux tube than when the quarks arc sur- 
rounded by gluons fixed to the Coulomb gauge or when 
they are connected by a thin string of gluons. Although 
in ref. [14] only mesonic states were considered, gener- 
alization to baryons is straightforward, with the gluonic 
strings attached to each of the three quarks contracted 
at a point. 

In this work we opt for the calculation of density- 
density correlators [15,16], 
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C(rMM)= / d'V (h\p u (v' +v,t 2 )p d (v',t 1 )\h) (7) 



where the density operator is given by the normal order 
product 



p u (r,t) = : w(r,t) 7o u(r,f) : 
so that disconnected graphs are excluded. 



(z,T/2) 



(8) 




(x n ,0) 



In the case of baryons three density insertions are 
needed and the correlator is given by 



C(n,r 2 ,t) = J d s r' (h\p a {v',t)p u {v' + r 1) t)p"(r' + v 2 ,t)\h) 

(10) 

where we have taken all the density insertions to be at 
equal times. This involves two relative distances. Like 
(7), it can be computed efficiently by Fast Fourier Trans- 
form (FFT). 





-r^ (x 0) (z,T/2) ^(xo,0) 



FIG. 1. Density-density correlator for a meson, ti, t2, 
T/2 — ti and T/2 — are taken large enough to isolate the 
mesonic groundstate. 



For mesons the density-density correlator is shown 
schematically in Fig. II. Inserting a complete set of 
hadronic states we find 



rii ,n f ,n p,q 



e -E nf (p)(T/2-t 2 ) 

K f (p) 

e -B n (p+q)(t 2 -t 1 ) 



(n f ,p\p |n,p + q)e 



(n, p + q\p \m,p) 



s«(p + q) 



E ni (p) 



■(m,p\h) ■ (9) 



The sum over all excitations of the initial, rii, and final, 
rif, states yields the ground state hadron in the limit 
(T/2 — t 2 ) — » oo and t\ — > oo, where T is the lattice ex- 
tent in Euclidean time. Since we take periodic boundary 
conditions the maximum time separation is T/2. Enforc- 
ing zero momentum would require a summation over the 
spatial volume on the source or sink site. This is techni- 
cally not feasible since it involves quark propagators from 
all to all spatial lattice sites. Instead, the suppression 
of the non-zero momenta and other higher excitations is 
obtained by choosing the largest possible time separation 
from the source and the sink. We take the same spatial 
coordinate for the source and the sink, namely xo = z. 
As a starting point we take, in this work, density inser- 
tions to be always at equal times. A disadvantage of 
the density-density correlators is that they are subject to 
more severe finite size effects, having typically larger spa- 
tial extent (~ twice) than Bethe-Salpeter amplitudes [6]. 




(z,T/2) 



(x ,0) 



FIG. 2. The upper two diagrams show the three-density 
correlator for a baryon as defined in Eq. 10. As in Fig. 1 
the density operators are inserted far enough from the source 
and the sink so that the baryonic groundstate of interest is 
isolated. The lower diagram shows the equivalent two-density 
correlator after integration of one of the relative distances. 

The relevant diagrams are shown in the upper part of 
Fig. 2 for the nucleon or the A + , for the general case of 
density insertions at three unequal times. In addition to 
one density insertion for each it-quark line, one can have 
two density insertions on the same it-quark line. Evalua- 
tion of this second diagram requires the quark propagator 
G(r 2 ,ri) for all partial distances of the two arguments. 
This is beyond our present resources, and this diagram 
is not included. To check that the first diagram that we 
calculate provides by itself a reasonable description of 
the baryon wave function, we also compute the baryon 
wave function with two density insertions for the it and 
d quarks as shown in the lower part of Fig. 2. This may 
be viewed as the square of the one-particle wave function 
obtained from the full wavefunction by integrating over 
one relative coordinate. If the contribution of the second 
diagram is small then we expect that 
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y'd 8 raC(ri,r a ,t)~ J d 3 r' (h\p d (r', t)p u (r' + n, t)\h) 

(11) 

will be satisfied, even when the l.h.s of the equation is cal- 
culated using only the diagram with one density insertion 
on each quark line. This comparison will be performed 
in Sec.IV, B. 



III. RELATIONS AMONG THE DEFORMATIONS 
IN DIFFERENT CHANNELS 

The interpolating field for a rho meson is taken to 
be J M (x) = d(x)j l _ l u(x). The physical states of spin 
projections and ±1 for the rho are obtained us- 
ing interpolating fields Jo(x) = d{x)^u{x) and J± — 
d(x) [(71 =F i7a)/2] u(x) respectively. Due to rotational 
invariance the correlator 



C ss (r) 



E 



(J s (z)/(x) p«(x + r)jt(0)) 



(12) 



satisfies the relation Cn(r.ei) = C22(r.e 2 ) = C33(r.e 3 ) 
where &j is a unit vector along the j-axis. This means 

that for the ±1 channels we have (C + C++)/2 = 

C\\ + Ci2- Due to parity symmetry C = C++ which 

is satisfied after ensemble averaging and leads to the re- 
lations 



C++ — C__ — Cn 



a 



22 



(13) 



If we denote by T(r) the tranverse and by L(r) the 
longitudinal projection of C3s(r) with respect to the spin 
axis, the deformation a(r) is defined by L(r)/T(r) = 
1 + a(r). The deformation in the ±1 channels is then 
2T(r)/(T(r) + L{r)) ~ 1 - a(r)/2 for small a i.e if the 
spin-0 state of the rho is elongated (prolate) along the 
z-axis then the spin ±1 states will be "fiat" (oblate) by 
approximately half this amount. This observation is con- 
sistent with the data shown in Fig. 3. 

Non-relativistically one can use the Wigner-Eckart the- 
orem to obtain a similar result. The deformation is ob- 
tained by measuring the quadrupolc moment defined by 

Q = (JM\2(z 2 -^±^)\JM) 

— {JM\2r 2 Y 2Q {e,4>)\JM) 



J—M 



J 23 
-MOM 



(J\2r 2 P 2 (cos(d))\J) (14) 



for a state \JM) where J is the spin of the state, M the 
spin projection along the z-axis and (Jmq J m) * s ^ e ^"J 
symbol. Therefore if we know the deformation for one 
spin projection we can relate it to the rest. Evaluating 
the 3-j symbol for the quantum numbers of the rho we 
find that the deformation for M = is twice and opposite 



in sign to that for M = ±1. We will thus only show 
results for the spin-0 state, which exhibits the largest 
deformation. 

1.1 . , , , , 

4 Jz=o 
0.9 
0.8 
0.7 




FIG. 3. Density-density correlator for the rho in a definite 
spin state J z = (top) and J z — +1 (bottom), for quark 
separations transverse (line with Z label) and longitudinal 
(line with X and Y labels) with respect to the spin axis. The 
hopping parameter is k = 0.154. 



The interpolating fields for projecting to the physical 
A + spin states were given in Eq. 4. Just like for the de- 
formation of the rho in different spin projections, similar 
relations can be obtained for the physical components of 
the A + . The interpolating fields for positive spin pro- 
jections 3/2 and 1/2 are no longer related to the nega- 
tive ones since they involve different spinor components. 
However the cross terms of the type 7172 contribute an 
order of magnitude less to the correlator and if we neglect 
these terms we find for the density-density correlator of 
the ±3/2 state that C 3/2 ~ C_ 3/2 ~ 1 - a(r)/2 * where 
a(r) is the deformation if we take 73 in Eq. 4. 

A relation among all the spin projections of the A + 
can be obtained in the non-relativistic limit by applying 
the Wigner-Eckart theorem Eq. 14. We find that the de- 
formation for the M = ±3/2 states is equal and opposite 
to that for the M — ±1/2. Among the physical states 



*The small difference between C3/2 and C_3/2 is due to our 
limited statistics. 
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we will therefore show results only for the +3/2. Since 
we expect the deformation to be maximal when a 73 is 
used in Eq. 4 for the interpolating field, we will also look 
for deformation in this channel in addition to the +3/2 
physical state. In the unqucnchcd data where we observe 
a small deformation these relations between the various 
amplitudes, as far as the relative signs are concerned, are 
indeed satisfied. 



IV. QUENCHED LATTICE RESULTS 
A. Density-density correlation functions 

We have analysed 220 quenched configurations at f3 = 
6.0 for a lattice of size 16 3 x 32 obtained from the NERSC 
archive [17], using the Wilson Dirac operator with hop- 
ping parameter k = 0.15, 0.153, 0.154 and 0.155. The ra- 
tio of the pion mass to the rho mass at these values of k is 
0.88, 0.84, 0.78 and 0.70 respectively. Using the relation 
2am q = 1/k — 1/k c , with the critical value k c — 0.1571, 
we obtain for the naive quark mass m q values of about 
300, 170, 130 and 90 MeV respectively, where we used 
a- 1 = 1.94 GeV (a = 0.103 fm) from the string ten- 
sion [18] to set the scale. Alternatively, the scale could 
be set from the rho mass in the chiral limit. This ap- 
proach yields a~ x =23 GeV (a = 0.087 fm) [19,20], with 
a systematic error of about 10% coming from the choice 
of fitting range and chiral extrapolation ansatz, which is 
about twice as large as the statistical one [19]. We note 
that although that these simulations used the lattices of 
increasing bigger volumes the scale was not affected. In 
our discussion of quenched data we will use the value of 
a determined from the string tension. However, to com- 
pare the quenched with the unquenched results we will 
use the value extracted from the rho mass in the chiral 
limit since this determination is applicable both in the 
quenched and in the unquenched theory. 

We fix the source and the sink for maximum separation 
at t{ = a and tt = 17a, where a is the lattice spacing. 
The density insertions are taken in the middle of the time 
interval i.e. at t — 9a. To check that the time interval 
\t — ti\ = \t — tf \ = 8a is sufficient, we have performed an 
analysis on 27 configurations at k = 0.153, varying the 
time t where the density operators are inserted. 

The results for the rho correlator are shown in Fig. 4(a) 
for four different insertion times at k = 0.153. As can be 
seen, density insertions at t = 8a and 9a give the same 
correlator, reassuring us that the time separation is large 
enough to isolate the groundstate of the rho. In the other 
channels the results are similar, with larger statistical 
noise for the baryons. In Fig. 4(b) we show in addition 
the correlator for the spin-0 state of the rho, for quark 
separations along the spin axis (z) and perpendicular to 
it (x). As expected, the deformation is the same when 
measured from density insertions at t — 8a and t = 9a, 
since the groundstate is isolated by then. More surpris- 



ingly, this z — x asymmetry is almost unchanged when 
measured at t = 4a, even though the z- and x-profiles 
change appreciably. These findings indicate that the de- 
formation which we observe in more detail below is a 
robust, physical property of the rho meson in its ground- 
state as well as its low-lying excited states. 

(a) (b) 




2 4 6 8 10 12 14 2 4 6 8 
r/a r/a 

FIG. 4. (a) The density-density correlator for the rho, mea- 
sured from density insertions at time 9a, 8a, 6a and 4a. To 
avoid cluttering, the data as a function of r are averaged over 
bins of size Ar = 0.07a. (b) The density- density correlator 
for the rho, measured from density insertions at time 9a, 8a 
and 4a, for quark separations along the z- and a;-axes. In 
both cases, 27 configurations are used, and statistical error 
bars are omitted for clarity. 

In Fig. 5 we collect the correlators for the pion, the 
rho, the nucleon and the A + for the four different quark 
masses (k values) considered. We confirm an observation 
made in earlier studies [3,4] that the wave functions are 
not very sensitive to the bare quark mass. This behaviour 
is as expected in the bag model, but is inconsistent with 
nonrelativistic quark models where a much stronger mass 
dependence is predicted. Insensitivity on the bare mass 
can be understood from the consideration that the quarks 
are dressed and their effective mass is therefore not very 
much affected by changing the bare mass. From Fig. 5 
we also see that the rho wavefunction depends more on 
the quark mass than the other three. In particular the 
nucleon and A + wave functions hardly change as we go 
from k — 0.154 to K = 0.155 which corresponds to re- 
ducing the naive quark mass from m q ~ 130 MeV to 
m q ~ 90 MeV. 
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FIG. 5. Density-density correlators, C(r), for (a) the pion, 
(b) the rho, (c) the nucleon and (d) the A + versus |r| at 
k = 0.15,0.153,0.154 and 0.155. Errors bars are omitted for 
clarity. 

A direct comparison of the sizes of the four hadrons can 
be made in Fig. 6(a) where we plot the density-density 
correlators for k — 0.154. The pion has the smallest size, 
approximately half that of the A, whereas the rho has a 
size comparable to the nucleon and the A. 
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FIG. 6. (a) Density-density correlators, C(r), for the pion, 
the rho, the nucleon and the A + at k = 0.154 vs |r|. (b) 
Same as (a) but with two dynamical quarks at k = 0.157. 
The dynamical results will be discussed in section V. We used 
the rho mass to set the scale in both the quenched and the 
unquenched theory. Errors bars are omitted for clarity. 

Any asymmetry with respect to the spin axis z is best 
seen by comparing the correlator C(r) for r = (x,y, 0) 
and (x,0, z), i.e. in a plane perpendicular to the spin 
axis and a plane containing it. The resulting contour 
plots are shown for n = 0.154 in Figs. 7 and 8 for the pion 
and the rho and in Figs. 9 and 10 for the nucleon and 
the A + . The cigar shape is clearly visible in the case of 
the rho, whereas the pion and the nucleon look spherical 
(up to small lattice distortions) as expected since their 
spectroscopic quadrupole is zero. 

A small asymmetry appears for the A + , but it is not 
statistically significant. This is true for all the A + spin 
projections. In Fig. 10 we select the unphysical state 
obtained using interpolating field J3 because it should 
show maximal deformation. In the case of the A, quench- 
ing removes only part of the pion cloud, contrary to the 
case of the rho where it is removed completely. Still, 
our negative finding does not rule out a deformation of 
the A induced by pions. It may be that the asymmetry 
only shows up when the quark mass is further decreased 
so that the pion has a mass close to its physical value, 
although we do not observe a statistically significant in- 
crease in the deformation as we go from quenched quarks 
of naive quark mass m q ~ 300 MeV to m q ~ 90 MeV. It 
could also be that the asymmetry is enhanced and only 
becomes visible in full QCD, with a complete pion cloud 
made of light enough pions. This issue is addressed in 
the next section. 
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FIG. 7. Contour plot of the pion correlator, C(r), at 
k = 0.153, when r lies in the xz-plane (solid lines) or in the 
xy-plane (dashed lines) of size 16 2 . 




FIG. 8. Same as Fig. 7 but for the rho for J z = 0. We have 
included a circle to guide the eye. 

A more quantitative determination of meson deforma- 
tions can be obtained by computing the second moments 
of the quark separation r along the three axes. Fig. 11 
shows the second moments along the spin-axis, (z 2 ), and 
the average of the moments along the two transverse axes, 
((x 2 +y 2 )/2), plotted versus the pion mass squared. Here 
we use the rho mass to convert to physical units in or- 
der to be able to compare with the corresponding un- 
equenched results. 
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FIG. 9. Same as Fig. 7 but for the nucleon. 




FIG. 10. Same as Fig. 8 but for the A + . The interpolating 
field is J3. 

The spin-0 state of the rho shows an elongation along 
the spin-axis, which increases as the quark mass is de- 
creased. As already mentioned, since this is a quenched 
calculation, this deformation is not due to the pion cloud. 
Therefore the situation may change if dynamical quarks 
are included. This is studied in Section V. 

From the second moments we can obtain the charge 
root mean square (rms) radius of the mesons defined in 
the quark model by 

<r 2 h > =]T eg ((r 9 -R) 2 ) 
9 
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E q e q f<Pr (r/2) 2 C(r) 
/ d 3 r C(r) 



(15) 



where R is the coordinate of the center of mass and e q is 
the electric charge of the quarks. In the chiral (quenched) 
limit we estimate for the pion ^fr 2 ~ 0.35 fm using the 
rho mass to set the scale and 0.42 fm using the string 
tension [18] to be compared with the experimental value 
of 0.53 fm [21]. For the rho, using the rho mass and the 



string tension we find respectively 



0.37 fm and 



0.44 fm. The ratio yr 2 /r 2 ~ 1.06 from the lattice data 

is somewhat small compared to the value 1.15, obtained 
from the experimental value 0.53 fm for the pion [21] and 
0.61 fm for the rho [22] calculated in a Dyson-Schwingcr 
equation approach using /„. = 93 MeV. This discrep- 
ancy between the experimental and lattice results may 
be due to two reasons: 1) we are using the quenched the- 
ory where the pion cloud is eliminated and 2) we are far 
from the chiral limit and making a linear extrapolation 
in m 2 may be problematic, especially since chiral loops 
give a logarithmically divergent contribution to the pion 
and proton charge radii [23]. Lighter quark masses will 
be required to check the chiral extrapolation. In section 
V we will examine the effects of unquenching. 
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FIG. 11. (z 2 ) and ((x 2 + y 2 )/2) versus the pion mass 
squared in physical units. The lines are linear fits to the 
quenched data. The pion (dash-dotted line) is spherical. It 
is smaller than the rho, whose transverse size with respect to 
the spin axis (solid line) is smaller than its longitudinal size 
(dashed line). The rho (in the spin-0 projection) is shown to 
be cigar shaped, particularly in the chiral limit. The circles 
show full QCD results, discussed in section V. 

A similar comparison of (z 2 ) with ((x 2 + y 2 )/2) for the 
A + , using as interpolating field J3 defined in Eq. 4, gives 



(z 2 ) > ((x 2 +y 2 )/2) (the reverse is true for the states with 
spin projection ±3/2), giving suspicion of a deformation. 
However, statistical errors are larger than the signal, so 
that we cannot claim to see any significant deformation. 



B. Three-density correlation functions 

We have analysed 30 configurations at k = 0.15 and 
k = 0.154. Since we now have to consider two relative dis- 
tances the calculation of the three-density correlations is 
more demanding even using FFT. As explained in section 
II we only compute the diagram with density insertions 
on different quark lines. However, we can check the qual- 
ity of this approximation by integrating our three-density 
correlation over one relative distance, and comparing the 
result with the two-density correlation studied in the last 
subsection. If the three-density correlation contains com- 
plete information, both expressions should agree as per 
Eq. 11. 

This comparison is performed in Fig. 12 for the nu- 
cleon, for two quark masses. Although the three-density 
correlation, as expected, is subject to larger finite-size ef- 
fects visible for large quark separation, both methods give 
virtually identical results, indicating that we captured 
the dominant contribution to the three-density correla- 
tion. A similar conclusion holds for the A. 

In any case, the usefulness of three-density correlators 
is to supplement two-density correlators and expose more 
detailed structure: single-particle observables, such as 
the quadrupole moment, can be extracted directly from 
the two-density correlators. 
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FIG. 12. Single particle density for the nucleon (a) for 
k = 0.15 and (b) for « = 0.154. 
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FIG. 13. u- and d- quark spatial distributions in the neu- 
tron for k — 0.15 and k = 0.154. The errors bars for the 
Tdu = case are comparable to those for the rdd = case and 
are omitted for clarity. 

The spatial u- and d- quark distributions in the nu- 
cleon can be investigated by fixing the relative position 
of the other two quarks. We note that, since we are using 
degenerate u— and d— quarks, only the electric charge 
differentiates the proton from the neutron. In Fig. 13 
we show the spatial distributions of the u and d quark in 
the neutron, when the relative distance between the other 
two quarks is fixed to zero. For both quark masses consid- 
ered, the d-quark spatial distribution is slightly broader 
than that of the u-quark. Since the total charge of the 
two d-quarks is -2/3 and that of the u-quark +2/3, the 
broader d-quark spatial distribution indicates that the 
charge root mean square radius of the neutron is nega- 
tive. For the A + the two distributions are the same. The 
charge radius squared can be evaluated using 



'ch/ 



fd^fd^Zl^e^jv^) C(n,r 2 ) 
J<P ri J<pT 2 C( ri ,r 2 ) 



(16) 



where r^(ri, r 2 ) is the distance of each quark to the cen- 
ter of mass, in terms of the relative distances ri and r 2 . 
Estimating this integral by a discrete lattice sum we ob- 
tain for the proton charge rms (r p /a) 2 = 20 ± 3.5 or 

- 0.45 ± 0.04 fm at « = 0.15 and (r p /a) 2 = 27 ±5.5 



hi' jrp ~ 0.52 ± 0.06 fm at k — 0.154. If we extrapo- 
late these two values linearly in m 2 to the chiral limit we 
find ~ 0.59(4) fm, compared to the experimental value 
0.81 fm. Again we have used a -1 = 1.94 GeV as de- 
termined from the string tension to convert to physical 
units. Using the nucleon mass in the chiral limit to set 
the scale gives a -1 = 1.88(7) GeV, very close to the value 
obtained from the string tension. If instead we would use 



the rho mass to set the scale, then a 1 w 2.3 GeV which 
gives for the proton, in the chiral limit, a smaller value 

equal to ^Jr 2 ~ 0.50(4) fm. It has been pointed out that 

chiral logs can affect the chiral extrapolation of the radii 
and increase their value [23] . To check whether the chiral 
logs will produce larger values closer to the experimental 
result, additional k values closer to the chiral limit will 
be needed. The neutron charge radius square r 2 comes 
out negative for both values of k as expected. We ob- 
tain r 2 n /r 2 p 0.25 ± 0.08 and -0.29 ± 0.12 at k = 0.15 

and k = 0.154 respectively. These values are consistent 
with the experimental value of —0.146 albeit with large 
errors. Computing as a check the same quantity for the 
A + , we obtain zero within error bars, in agreement with 
our earlier observation that the d- and it- quark spatial 
distributions in the A + are the same. 
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FIG. 14. d-quark spatial distribution with respect to the 
center of mass of the two u-quarks, for different u — u separa- 
tions (a) in the nucleon spin +1/2 projection and (b) in the 
A + spin +3/2 projection (k = 0.15). The straight lines mark 
the value of the wave function when r uu — and the u — d 
quark separation is 0.2, 0.4 and 0.8 fm, to be compared with 
T'cm = and r uu = 0.2, 0.4 and 0.8 fm respectively. 

Instead of fixing the relative distance between the two 
w-quarks in the proton or A + to zero, we can fix it to 
various non-zero values. In Fig. 14 we fix the u — u sep- 
aration, r uu , (between and 8a along a principal lattice 
axis) , and show the distribution of distances r cm between 
the d-quark and the u — u center of mass. Detailed in- 
formation about nucleon structure is contained in this 
figure. 

In particular, let us compare configurations (uu) — d 
where the two u-quarks at at the same place i.e. r uu = 0, 
and u — d— u where the d-quark lies in the middle of the 
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two u-quarks, i.e. r cm = where r cm is the distance from 
the d-quark to the center of mass of the two it-quarks. 
As we will see in the next paragraph, effective models 
(both Y- and A-Ansatze) predict equality of the wave 
functions among these two configurations provided they 
have the same total size (ie. r cm in {uu) — d equal to r uu 
in u — d — u). Inequality reveals a finer structure than 
captured by effective models. It is clear from Fig. 14 
that the top horizontal line {{uu) — d with r cm = 0.2 
fm) falls below the + sign at r cm = (u — d — u with 
r uu — 0.2 fm), indicating a relative suppression of the 
{uu) — d configuration, i.e. a mutual u — u repulsion. The 
same effect is visible for a larger configuration of size 0.4 
fm (second horizontal line), but goes away or even gets 
inverted for yet larger configurations. Somewhat surpris- 
ingly, this u — u repulsion at close range seems weaker 
in the A + where the two u-quarks are in the same spin 
state (Fig. 14(b)) than in the nucleon where they can be 
in opposite spin states (Fig. 14(a)). This example, where 
a quantitative refinement could be obtained straightfor- 
wardly by considering lighter quarks on a larger lattice, 
serves to illustrate the wealth of information contained 
in three-density correlators. 

Having, in the non-relativistic limit, the wave func- 
tion of a baryon it is interesting to ask whether we can 
deduce the potential which would yield this wavefunc- 
tion. The relevant potential will of course depend on the 
two relative distances between the quarks. The issue is 
whether we can further reduce the degrees of freedom to 
effectively write the confining potential in terms of one 
distance only. Two such proposals exist in the literature, 
which make different predictions for the linear rise of the 
baryonic potential, (z) The so called Y-Ansatz can be 
derived by a strong coupling argument [24]: the baryon 
potential grows in proportion to the minimal length of 
gluonic string necessary to join together the three static 
quarks. The three strings join at the Steiner point, {ii) 
The so called A-Ansatz is derived from a center vortex 
picture of confinement [11]. The baryonic potential is 
simply a sum of two-body qq potentials and grows like 
§( r i2 + ?"23 + f"3i)) where a is the qq string tension and 
Tij the qi — qj distance. The difference between these 
two Ansatze depends on the geometry of the three quark 
system. For instance it vanishes when the three quarks 
are aligned, since then the Steiner point which minimizes 
the length of flux tube joining the three quarks coincides 
with the position of the middle quark. Unfortunately this 
difference never exceeds ~ 15%, so that very accurate re- 
sults at large quark separations ~ 1 fcrmi are needed to 
ascertain which model, if any, is correct. In a recent lat- 
tice calculation of the baryonic potential it was shown 
that the A-Ansatz is favoured for distances up to 0.8 fm 
[25]. Although this finding is in agreement with the re- 
sults of [26], a different analysis of similar lattice results 
has led others [27] to the conclusion that the Y-Ansatz 
is preferred. This issue is currently under further study 
both on the lattice [28] and phenomenologically [29] . 
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FIG. 15. Nucleon wave function (a) versus half the 
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FIG. 16. Theta dependence of the nucleon wave function 
when the two relative distances are fixed to 0.6 ± 0.001 fm. 
The filled squares are the central values of the wave function. 
The solid line is the fit using the A-Ansatz for the angular 
dependence and the dashed line is the corresponding fit using 
the Y— Ansatz. Note that for cos (9) < —0.5 the Y-Ansatz 
coincides with the A-Ansatz. 
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To examine whether our wave function results favour 
one of these two Ansatze, we plot them as a function of ry 
(Fig. 15(a)) and of r& (Fig. 15(b)), where ry is the min- 
imal total length of the flux strings and ta = l/2(ri2 + 
r 23 +F31). The scatter of the data is visibly smaller when 
plotted as a function of ta , indicating that r& is a better 
effective variable than ry. Furthermore, consider the so- 
lution to the non-relativistic Schrodinger equation with 
potential Vy oc ry or Va k ta- It is an Airy function 
which asympotically decays as exp(— cr 3 / 2 ). Fits of the 
wave function of the nucleoli to this asymptotic form are 
shown in Fig. 15. The fit using the r& as the relevant 
distance yields a x 2 /d.o.f = 0.4 whereas using ry one 
gets % 2 /d.o.f = 1.0. Therefore both Ansatze provide a 
surprisingly good description of the wave function, with 
a preference for the A-Ansatz. 

As a more stringent test, we fix the relative distances 
I rx I = I r2 1 between the d and the two u quarks, and 
study the wave function as a function of the udu an- 
gle, cos(6*) = n.r2/|ri||r2|, that the relative distances 
make. Both for the A- and Y— Ansatze the wave func- 
tion should only depend on this angle but with a different 
functional dependence. In Fig. 16 we show the wave func- 
tion for the nucleon versus the cosine of the angle. Again 
the A-Ansatz provides a better description of the data. 



V. DENSITY-DENSITY CORRELATORS WITH 
TWO DYNAMICAL QUARKS 

As we noted in the previous section, if the pion cloud 
is responsible for hadron deformation, then the quenched 
and unquenched results may differ significantly. To in- 
vestigate the importance of dynamical quarks, we used 
the SESAM [13] configurations with two dynamical de- 
generate quark species at /3 = 5.6 on a lattice of size 
16 3 x 32. The lattice spacing determined from the rho 
mass in the chiral limit is a -1 = 2.3 GeV [13] which is 
same as for the quenched theory at (3 = 6.0 and therefore 
the physical volume is the same as in our quenched calcu- 
lation. We use this determination of the lattice spacing, 
which is applicable in both the quenched and the un- 
quenched theory, to compare results among both. We 
have analysed 150 configurations at K = 0.156 and 200 
at k = 0.157. The ratio of the pion mass to rho mass is 
0.83 at k — 0.156 and 0.76 at k = 0.157. These values are 
close to the quenched mass ratios measured at k = 0.153 
(0.84) and re = 0.154 (0.78) respectively, allowing us to 
make pairwise quenched-unquenched comparisons. 

In Fig. 6(b) we plot the two-density correlators for the 
four hadrons for light sea quarks (k = 0.157). Com- 
paring with the quenched results (re — 0.154), one sees 
that the pion size remains the same, but that the rho, 
nucleon and A + sizes increase. The A + is now clearly 
larger than the rho and the rho correlator decays more 
slowly than the nucleon. This slower decay is also visible 
in the quenched case for k = 0.154 and 0.155 but not for 



the smaller values of re. The invariance of the pion size 
can be also seen in Fig. 11 which shows the second mo- 
ments. This result is consistent with calculations in the 
Dyson- Schwinger framework which concluded that the 
pion cloud, consisting of physical pions, contributes only 
15% to the root mean square radius of the pion [30]. In 
our simulations the effect is expected to be even smaller 
since the pion mass is larger than 600 MeV. On the other 
hand, Fig. 11 shows that unquenching increases the rho 
size and therefore the ratio of the rho charge radius to 
the pion's approaches the experimental value. This ra- 
tio is expected to be more reliable than the individual 
moments since lattice artifacts partly cancel out. 
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FIG. 17. Asymmetry for the rho wave function (a) in the 
quenched approximation at re — 0.153 and (c) at re = 0.154 
and (b) for two dynamical quarks at re = 0.156 and (d) 
re = 0.157 . 

From Fig. 17 we compare the quenched and un- 
quenched asymmetry in the rho for comparable ratios 
m v /m p . The main observation is that the asymmetry in 
the rho grows in full QCD. This is clearly seen in the (a)- 
(b) comparison, i.e. for the heavier dynamical quarks. 
For the lighter quarks (c)-(d), this growth in the asym- 
metry is still there but the effect seems less pronounced. 
On the other hand, finite-size effects are more important. 
Disentangling the two would require larger lattices. 
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FIG. 18. Three-dimensional contour plot of the correlator 
(black): upper for the rho state with spin projection (cigar 
shape) and lower for the A + state with +3/2 (slightly oblate) 
spin projection for two dynamical quarks at re = 0.156. Values 
of the correlator (0.5 for the rho, 0.8 for the A + ) were chosen 
to show large distances but avoid finite-size effects. We have 
included for comparison the contour of a sphere (grey). 

The same analysis for the A + gives no definite results 
regarding the deformation. Whereas in all channels we 
do observe a deformation obeying the sign relations ob- 
tained in Section II with the Wigner-Eckart theorem, the 
deformation remains small with large statistical errors. 
However it is interesting to look at a three-dimensional 
contour plot for the A + in the +3/2 spin state and com- 
pare it to the rho 0-spin state. In Fig. 18 we show these 
contour plots for the heavier sea quarks that we analyzed. 
The elongation in the rho is clearly visible whereas the 
A + appears to be squeezed. Note that squeezing in this 
channel means that the unphysical channel using an in- 



terpolating field a J3 in Eq. 4 is elongated. The statistical 
uncertainties are larger for the lighter sea quark but the 
trend for the deformation remains the same. The overall 
conclusion is that unquenching seems to increase the de- 
formation for the p and the A, which would imply that 
the pion cloud contributes to the deformation of these 
hadrons. However more statistics, a bigger lattice and 
lighter quark masses will be needed to consolidate this 
observation. 



VI. CONCLUSIONS 

Two-density correlators for mesons and baryons are 
shown to contain rich information on their structure. For 
the baryons these correlators reduce to the one-particle 
density in the non-relativistic limit. 

The overall quark mass dependence of the correlators 
for the pion, the rho, the nucleon and the A + is found to 
be rather weak. The rho correlator shows the strongest 
dependence on the quark mass, and the nucleon and the 
A the weakest. For the quark masses that we have stud- 
ied in this work unquenching has the strongest effect on 
hadron sizes for the rho and the A, and the weakest for 
the pion and the nucleon. In the quenched approxima- 
tion we find that the rho state with 0-spin projection is 
prolate, whereas for the A + no statistically significant 
deformation is seen. Since for the rho there are no vir- 
tual pions from backward moving quarks in the quenched 
approximation, the rho deformation is a dynamic prop- 
erty of quarks forming a vector particle. Adding sea 
quarks while keeping the pion to rho mass ratio con- 
stant, we observed an increase in the rho deformation, 
and a slight oblate deformation of the A + +3/2 spin 
state. The pion cloud may be responsible for this ef- 
fect. As already pointed out, lighter sea quark masses 
together with a larger lattice to well contain the hadrons 
will be needed in order to study further the role of the 
pion cloud. 

Three-density correlators for baryons, although com- 
puted neglecting insertions on the same quark line, are 
shown to reproduce the two-density correlators when in- 
tegrated over one relative coordinate. Therefore one ex- 
pects the contribution of the neglected diagram with two 
insertions on the same quark line to be small. With the 
three-density correlators baryonic structure can be ex- 
plored in more detail. In the neutron we clearly detect a 
broader <i-quark spatial distribution compared to that of 
the it-quark. This accounts for the negative charge square 
radius of the neutron observed experimentally. By com- 
paring the u- and d- spatial distributions in the proton 
and in the A + we observe that there is a preference for 
the two it-quarks to be at 180° rather than at the same 
place. More statistics are needed to consolidate the trend 
observed here. 

Information on the baryonic potential can be extracted 
from fits to the three-density correlators of the nucleon. 
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By performing a fit to the radial and angular dependence 
we find that the baryonic wave functions are better de- 
scribed by a confining potential which is the sum of two- 
body potentials known as the A-Ansatz, at least for rel- 
ative distances of <~ 1 fm that we can probe in this work. 
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